---
title: "Van Egeren et al (2018) Figs. 1-2"
output: html_notebook
---

```{r include=FALSE}
library(readr)
library(ggplot2)
library(reshape2)
library(RColorBrewer)
source("./utils.R")
```

## Fig. 1b

```{r include=FALSE}
filepath_dists = "./dists/"
filename = paste(filepath_dists,"cd8_norm.txt", sep="")
temp = read_tsv(filename, skip=1, col_names=FALSE)
cd8_dist = melt(temp)
filename = paste(filepath_dists,"l1210_norm.txt", sep="")
temp = read_tsv(filename, skip=1, col_names=FALSE)
l1210_dist = melt(temp)
```
```{r echo=FALSE}
ggplot() + geom_line(aes(x=800:1200/(1000), y=dlnorm(800:1200/(1000), meanlog=mu, sdlog = sqrt(sig))), color="orange", size=0.8) + xlim(0.85, 1.15) + geom_line(aes(x=cd8_dist$value * sqrt(sig2) + 1), stat="density", color="purple", size=0.8) + theme_bw()+ geom_line(aes(x=l1210_dist$value * sqrt(sig2) + 1), stat="density", color="green3", size=0.8)+ theme(legend.position = "none", axis.title.x=element_blank(), axis.title.y=element_blank(), text = element_text(size=20), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"), panel.border = element_blank())
```

## Fig. 2b

```{r include=FALSE}
filepath = "./"
mean_files = c()
run = c()
lifetimes = c()
variance = c()
for (j in 1:10){
  filename = paste(filepath, "lognorm_mean_fit/heritable/mean_fit_sim_",j,"type_0.oevo", sep="")
  mean_files = c(mean_files, filename)
  
  variance = c(variance, -1)
  lifetimes = c(lifetimes, -1)
  run = c(run, j)

  for (i in c(1,5,10,50,100)){
    for (k in 4:4){
    run = c(run, j)
    variance = c(variance, k)
    lifetimes = c(lifetimes, i)
    filename = paste(filepath, "lognorm_mean_fit/", i, "/", k, "/mean_fit_sim_",j,"type_0.oevo", sep="")
    mean_files = c(mean_files, filename)
    }
  }
}

mean_data = do.call(rbind,lapply(mean_files, function(x){read_csv(x, skip=1, col_names=FALSE)}))

colnames(mean_data) = c("time", "mean_fit")
mean_data["variances"] = rep(variance, each=365001)
mean_data["lifetimes"] = rep(lifetimes, each=365001)
mean_data = mean_data[mean_data$time < 10000,]
```

```{r include=FALSE}
mean_data = mean_data[mean_data$time %% 100 == 0,]
mean_data["mean_fit_adj"] = (mean_data$mean_fit)
aggregated_traj = aggregate(mean_data$mean_fit_adj, by=list(mean_data$lifetimes, mean_data$variances, mean_data$time), FUN=mean)
colnames(aggregated_traj) = c("lifetimes", "variance", "time", "mean_fit")
aggregated_traj["sd"] = aggregate(mean_data$mean_fit_adj, by=list(mean_data$lifetimes, mean_data$variances, mean_data$time), FUN=sd)$x
aggregated_traj["colors"] = as.numeric(sapply(aggregated_traj$lifetimes, FUN=to_colorscale3))
```

```{r echo=FALSE}
min_y = 0.97
max_y = 1.5
#min_y = -Inf
#max_y = Inf
ggplot() + geom_line(data=aggregated_traj[aggregated_traj$lifetimes > 0,], aes(x=time, y=mean_fit, group=lifetimes, color=as.numeric(colors))) + geom_ribbon(data=aggregated_traj[aggregated_traj$lifetimes > 0,], aes(x=time, ymin=pmax(min_y, mean_fit-sd), ymax=pmin(mean_fit+sd, max_y), group=lifetimes, fill=as.numeric(colors)), alpha=0.3) + theme_bw() + scale_color_gradient2(low="#cccc00", mid="#ff3399", high="#9900ff", midpoint=4) + theme(legend.position = "none", axis.title.x=element_blank(), axis.title.y=element_blank(), text = element_text(size=20), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"), panel.border = element_blank()) + scale_fill_gradient2(low="#cccc00", mid="#ff3399", high="#9900ff", midpoint=4) + geom_line(data=aggregated_traj[aggregated_traj$lifetimes < 0,], aes(x=time, y=mean_fit), color="#4400cc") + geom_ribbon(data=aggregated_traj[aggregated_traj$lifetimes < 0,], aes(x=time, ymin=pmax(min_y, mean_fit-sd), ymax=pmin(mean_fit+sd, max_y)), alpha=0.3, fill="#4400cc") + xlim(0, 10010)
```

## Fig. 2c

```{r include=FALSE}
filepath = "./"
all_dat_twopt = c()
lifetimes = c()
input_vars = c()
for (k in 4:4){
  for (i in c(1,3,5,10,50,"heritable")){
    lifetimes = c(lifetimes, i)
    input_vars = c(input_vars, k)
    for (j in 1:1000){
      filename = paste(filepath, "lognorm_ss_dist/", i, "/", k, "/fit_sim_",j,"type_0.oevo", sep="")
      temp = read_csv(filename, skip=1, col_names=FALSE)
      est_mean = mean(unlist(t(temp[1,2:101])))
      to_deliver = unlist(t(temp[1,2:101]))/est_mean
      all_dat_twopt = c(all_dat_twopt, to_deliver)
    }
  }
}
dist_data = data.frame(fit=all_dat_twopt, lifetimes=rep(lifetimes,each=100000), input_variances=rep(input_vars,each=1000))
```

```{r echo=FALSE}
g1 = plot_dists(dist_data)
g1
```

## Fig. 2d

```{r include=FALSE}
filepath_dists = "./dists/"
filename = paste(filepath_dists,"cd8_norm.txt", sep="")
temp = read_tsv(filename, skip=1, col_names=FALSE)
cd8_dist = melt(temp)
filename = paste(filepath_dists,"l1210_norm.txt", sep="")
temp = read_tsv(filename, skip=1, col_names=FALSE)
l1210_dist = melt(temp)
all_dat_lgnorm_hold = dist_data[dist_data$lifetimes==10,]$fit

est_mean = mean(all_dat_lgnorm_hold)
to_plot = all_dat_lgnorm_hold/est_mean
#est_var = var(log(to_plot))
mu = mean(log(to_plot))
sig = mean(((log(to_plot) - mu)^2))
mu2 = mean(to_plot)
sig2 = var(to_plot)
```
```{r echo=FALSE}
temp = density(all_dat_lgnorm_hold, n=40)
ggplot() + geom_line(aes(x=800:1200/(1000), y=dlnorm(800:1200/(1000), meanlog=mu, sdlog = sqrt(sig))), color="orange", size=0.8) + xlim(0.85, 1.15) + geom_line(aes(x=cd8_dist$value * sqrt(sig2) + 1), stat="density", color="purple", size=0.8) + theme_bw()+ geom_line(aes(x=l1210_dist$value * sqrt(sig2) + 1), stat="density", color="green3", size=0.8)+ theme(legend.position = "none", axis.title.x=element_blank(), axis.title.y=element_blank(), text = element_text(size=20), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"), panel.border = element_blank()) + geom_line(aes(x=all_dat_lgnorm_hold), stat="density", linetype="dashed") + geom_point(aes(x=temp$x, y=temp$y), size=2, shape=15)
```
